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1. Introduction 


During the second year of the research program efforts have continued to concen- 
trate on the development of the numerical methods that will form the computational 
part of the turbulence closure scheme. Studies have continued on the wave model for the 
two-dimensional shear layer. This configuration is being used as a test case for the clo- 
sure schemes. Several numerical schemes for the solution of the non-separable Rayleigh 
equation have been developed. This solution is required for the closure scheme in more 
complex geometries. The most efficient method found is a hybrid scheme that combines 
both pseudospectral and finite difference techniques. In addition, conformal transforma- 
tion techniques have been developed to transform the arbitrary geometry of the jet to a 
simple computational domain. We have also continued our study of the shock structure in 
arbitrary geometry jets and multiple jets. These developments are described briefly below. 

2. Numerical Methods for Arbitrary Geometry Jets 

In this section we describe the progress we have made in the solution of the non- 
separable form of the Rayleigh equation. This efficient solution of this equation for ar- 
bitrary geometries is a key element in our turbulence closure scheme. The two areas in 
which progress has been made are (i) the development of a algorithm for the solution of 
the equation and (ii) the development of conformal transformation methods to transform 
the jet’s arbitrary geometry into a simpler computational domain. 

2.1 Numerical Algorithm: Pseudospectral Hybrid Scheme. 

A method has been developed to determine the eigensolutions of the Rayleigh equation 
in flows of arbitrary geometry. The equation to be solved is: 

(A-/? 2 )p+(2/?/n)VWVp = 0, (2.1) 
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with boundary conditions: 


p is finite and p — > 0 at infinity, 


( 2 . 2 ) 


where p is the pressure fluctuation, W(x,y ) is the axial mean velocity, (3 is the axial 
wavenumber, fl = u> — (3W, and oj is the wave frequency. The technique that has been 
developed to solve (2.1) and (2.2) is a “hybrid” scheme. The method consists of combining 
a finite difference approximation with a pseudospectral approximation. Consider a solution 
of the form: 

PN = 5^a(*y.y)/y(*), (2-3) 


where the summation is taken over 0 to N. Here the /y’s are functions defined by: 

(_!)<+■(! - x i)T'„(x) 


fi( x ) 


cjN 2 (x — Xj ) 


(2.4) 


where, 


Xj = cos (irj/N) for 0 < j < N, 


(2.5) 


and where lyy is the JV-th order Chebyshev polynomial. In addition, co = = 2 and 

Cj = 1 for 1 < j < N — 1. 


The discretization of (2.1) is obtained by substituting (2.3) and noticing that (2.4) 
and (2.5) imply: 

*<*>-{£ lit < 2 - 6) 

The resulting discretization is given by: 


{« dyy+H{x k ,y)dy}a(x k ,y ) - (3 2 a{x k ,y ) 

+ [{d xx + G{x k ,y)d x }^2a{x j ,y)f j {x)] X — Xk “ 0, 


(2.7) 


for 1 < k < N — 1. In (2.7) the functions H and G are defined ii terms of 2/?/f2 and t 1 
coiiipoiients of the gradient of W . 
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To outline the determination of the eigenvalues (3 for the Rayleigh problem consider 
the model shear domain: 

D:={{x,y)e R 2 - 1 < x , y < 1}. (2.8) 

Now, to determine the boundary conditions for the computational domain, D notice that 
the lines y ± 1 correspond to the edges of the shear layer. On the boundaries of the shear 
layer the axial mean velocity W becomes constant. Whence, (2.1) reduces to the Helmholtz 
equation that has a known analytic solution, say: 

p{x, ±1) = Rn{x, ±1), (2.9) 

where R n would involve Bessel or Hankel functions and exponentials in polar coordinates. 
In addition, for the lines x = ±1 a boundary condition based on the periodicity of the 
pressure fluctuation is assumed. For example of the model domain D , set p(±l,y) = 0. 
This could correspond to pressure fluctuations that are odd about both the major and 
minor axes of the jet. The determination of the eigenvalues proceeds as follows: 

i) The derivatives in (2.7) for the spectral direction are determined as functions of the 
collocation points, Xj. 

ii) Equation (2.7) is recast as a two-dimensional vector allowing the use of an explicit 
integration scheme in the y- direction. 

iii) The function Ro(xj, +1), from (2.9) and its derivative in the y-direction are evaluated 
at the collocation points on the line y = +1. This gives the initial values for the 
discretized equations that are then integrated to y = 0. 

iv) Step iii) is repeated for the functions i? n (xy,+l) for n < N — 2. 

v) Steps iii) and iv) are repeated, starting the integration at y = — 1 and using the 
functions R n (xj,— 1) as the starting conditions. 

vi) The integrated solutions are then matched at y = 0. The matching of the so’ at: is 
and their derivatives gives a 2 N X 2 N matrix, M(/?). 
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vii) The value of (3 satisfying det[M(/?)] = 0 for a fixed frequency is then computed 
by providing a first guess for (3 and then iteratively repeating steps iii) through vi) 
updating the value (3 after each iteration. 

A code implementing this “hybrid” scheme has been developed for 
circular and elliptic geometries. The code has been validated for the circular case 
and in some elliptic ases. Further verification of the elliptic jet case is depends on the 
coordinate transformations described below. For the circular jet a velocity profile of the 
form: 

iy(r) = (1 + tanh[(l - r)/20])/2, (2.10) 

has been assumed, where 6 is the momentum thickness of the mixing layer. The elliptic 
test case assumes an ellipse of aspect ratio 2 at the half velocity point of the mean velocity 
profile. 

2.2 Computational Domains for Jets of Arbitrary Geometry 

2.2.1. Introduction - An efficient method establishing the computational domain for 
jets of arbitrary geometry is currently being developed. The technique to be used involves 
generating conformal mappings which carry standard computational domains into the cross 
sections of a given jet. 

To begin the discussion, the general types of geometries encountered in jet flows will 
be outlined. Two topologically distinct cross sections occur in a jet. The first type of cross 
section corresponds to the annular shear region surrounding the jet’s potential core. Since 
the velocity in the potential core is approximately uniform, the computation needs only to 
be performed in the annular shear region. Therefore, annular domains constitute the first 
type of geometry in jet flows: doubly connected regions. 

The second general type of geometry encountered in jet flows corresponds to the region 
ft a. jet downstream of the pou ntial core. In this case, the computation must be performed 
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in the entire cross section of the flow. Thus, the second class of geometry are the simply 
connected, disc-like regions. 

The goal of the research on jet geometry is to establish coordinate transformations 
which map standard computational regions into cross sections of a arbitrary jet. The 
s'andard regions are chosen to support the numerical solution of the equations modeling 
jet flows. Since an arbitrary jet has two distinct types of geometries, two standard types 
of conformal maps and computational regions result. The following sections will outline 
the two types of conformal maps to be used and give an example for the doubly connected 
region of an elliptic jet. 

2.2.2. Simply Connected Regions - The standard computational region for this case is 
the unit disc centered at the origin of the complex plane. A conformal transform, F, map- 
ping the unit disc into a simply connected jet cross section will be generated numerically 
using a generalization of the method of Theodorsen ref. 1. To understand the geometry, 
consider Figure 1. Theodorsen ’s map carries the circle in the computational domain onto 
the outer edge of the shear layer in physical space. 

Let the boundary of the jet cross section be described in polar coordinates in physical 
space as: 

w — w{(j)) = R{4>) exp(i<f>), (2.11) 

where, R{4>) is assumed to be a 2n periodic function, and the region bounded by R(<j>) is 
required to be starlike with respect to the origin. From the development of Henrici ref. 2, 
Theodorsen’s method goes as follows: 

i) Assume the conformal map takes the form of: 

F(z) = zexp[h{z)], (2.12) 

where H(z) is an unknown function to be approximated. 
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ii) Evaluate (2.12) on the unit circle to obtain: 


3i//[exp(t0)] = ln{i2[<£(0)]}, 
ajj[exp(*0)] = <f>{0) - 0. 


(2.13) 


iii) Next, solve for the relation between the angled, in physical space and the angle, 0 , in 
computational space which may be shown to be expressed as: 


— 0 — K ln{i?[<£(0)]}, (2.14) 

where K denotes the principal value integral: 

k <f>{x) = -^ J q cot [( x - y)l 2 \ d y- (2.15) 

iv) Once equation (2.14) is known, equation (2.12) is used to generate the mapping from 
the disc to the jet cross section. 

Equation (2.14) is solved numerically by introducing the iterative equation: 

Afc+i(<^) = K ln{i2[0 + Afc(0)]}, (2.16) 

for A: = 0, 1, 2, ..., where A(0)^V(0) — 0. At each step an FFT is used to approximate 
(2.16). An initial guess for A k is given by 

A o {0) = 0. 


2.2.3. Doubly Connected Regions - The standard doubly connected computational 
domain is the circular annulus. The annulus will be normalized such that the outer bound- 
ary has radius one. The inner boundary radius, denoted by n, cannot be chosen arbitrarily 
but must be determined by the method used to construct the conformal map from the an- 
nulus to the jet cross section. The method developed by Theodorsen for simply connected 
domains was extended by Garrick ref. 3 to generate conformal transformations for annular 
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regions. The Theodorsen - Garrick scheme is the method that will be employed for doubly 
connected jet cross sections. The geometry is shown in Figure 2 

The Theodorsen - Garrick map, T, carries the outer boundary of the annulus onto the 
outer boundary of the jet cross section. Likewise, the inner annular boundary is assigned 
to the jet’s inner cross section boundary. 

To develop the map T, let 


Zi = exp(i<f>) for i = 0,1, (2-17) 

denote 2n periodic functions describing the boundaries of the jet cross section. It is required 
that the region, bounded by these functions, be starlike. Referring to Henrici ref. 2, the 
construction of T goes as follows: 

i) Assume that there exists a conformal map of the form: 

T(z) = zexp[G(z)], (2.18) 

where G{z) is the unknown function to be approximated, 
ii) Since G is defined on an annulus, it will have a Laurent series of the form: 

OO 

G{z) — a 0 ,o + 2 £ A n z n , (2.19) 

infty 

X^O 

where 

An = (ao,n ~ H n ai,n)/(1 ~ /* 2n ). (2-20) 

The a ti n’s, i = 0, 1, depend on the Fourier coefficients of the logarithm of the boundary 
functions given by (2.17). However, these Fourier coefficients depend on the angle, 0, 
in the computational domain, and may not be directly computed until a relation of 
the form <p{0) = f(0) is known. 
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iii) Therefore, a relationship between <f> and 6 must be determined. This is given by the 
pair of integral equations: 

MO) ~0 = H m ln{i2 0 [<£o (*)]} - G m ln^^ (*)]}, 

MO) - 0 = G m ln{iE o [^o(0)]} - H m ln^f^)]}. 

Here, 


H„ = K<£(0)+K ^(0), 

Q^{0) = K p4>(9), 

where K is defined as in equation (2.15) and 


k r m = - [* 

* J-1 r 


( 2 . 21 ) 


( 2 . 22 ) 


and where, 

OO 

kn(x) = ^(2// 2n /(l — /i 2n ) sin(nx). 
l 

Finally, n is given by: 

ix = exp ln{Ri[<f>i(0)]}/Ro{<j>o(0)]}do\ . (2.23) 

iv) Once equations (2.21) to (2.23) have been solved for <f>i and for fx the Fourier coef- 
ficients for ln{i? t [<^t(0)]} can be determined and substituted into equation (2.19) to 
give the Laurent series for G(z). Combining G with equation (2.18) then gives the 
desired conformal map. 


This set of nonlinear integral equations can be rapidly solved for the angle variations, 
<f>i, and for /x by using the iterative approximation to equations (2.21): 


^o,m+i 0 — H m ln[i?o(^o,m)] G m In[.Ri(</>i )m )], 

^ljm+i 0 — Gfj, ln[^Eo (<^o,m ) ] — ln[i?i (</>i, m )], 


(2.24) 


where for an approximated //. the operators defined by equations (2.22) are determined. The 
iterative scheme is very efficient because, all of the functions considered are approximated 
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using an FFT. At each step in the iteration a new value of p. is given by the exponential of 
the difference of the zero coefficients in the Fourier series, of the functions, ln{i?,[<£,, m (0)]}. 

2.2.4. General Remarks - The iterative schemes described in the sections above con- 
verge very rapidly for nearly circular regions. However, as the physical domains become 
more eccentric, i.e., for high aspect ratios, the methods converge slowly, and in some cases, 
not at all. This problem is fairly easy to overcome by using preliminary conformal trans- 
formations which carry the physical domain into a region which has a lower asp ct ratio. 
Moreover, all the iteration schemes may be recast using underrelaxation. Details of the 
underrelaxation technique for simply connected regions are described in Gutknecht ref. 4. 

It should be recalled that conformal maps locally preserve the angles between mapped 
curves. Therefore, an orthogonal coordinate grid in the computational domain will trans- 
late into an orthogonal grid in the physical space. Since the coordinates are orthogonal, 
the metric tensor associated with conformal maps will have only diagonal terms which are 
nontrivial. Therefore, conformal maps introduce a minimum number of derivative terms 
in the transformed equations of motion. 

2.2.5. Example - In this section an example of the conformal transform mapping the 
annulus onto a shear layer bounded by two confocal ellipses will be examined. The ellipses 
are defined by: 


x 1 + (4/3) j/ 2 = 1, 

(2.25) 

(16/9)x 2 + (16/5)y 2 = 1. 

(2.26) 


Equation (2.25) defines the outer edge of the shear layer, while equation (2.26) defines the 
inner edge. In order to estimate the accuracy of the map, an error measure will be defined 
for the points mapped to the boundaries of the cross section. Therefore, set 

Error Measure d = f |F[r(f£,0)] — 1 j , 
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where F denotes either equation (2.25) or (2.26) and R = 1 or fx depending on F. 
Let N denote the number of data points used to compute the FFT of the boundary 
functions, and let 6 denote the angle in the computational space. Calculations of 
the error measure are given in Table 1. 

The following figures, 3a and 3b, show a computational grid and its image 
in the physical domain for the numerical example exhibited here. For the figures: 
N = 1024 and the error is « O(10 -5 ). 

2.3 Shock Structure in Arbitrary Geometry Jets 

Models have also been developed to describe the shock structure and instabiltiy 
waves in single and multiple jets of arbitrary geometry. The initial results of this 
study were given by Morris, Bhat and Chen ref. 5. This paper, which is attached 
as an appendix to this report, described the use of a boundary element method to 
determine the shock spacing in jets of arbitrary geometry represented by a vortex 
sheet. This topic has been submitted for publication in the Journal of Sound 
and Vibration. At present we are extending the analysis to include the effects 
of finite mixing region thickness and the effects of the small-scale turbulence for 
both single and multiple jets. The results for the finite mixing layer thickness are 
described here. 

2.3.1 Effects of finite mixing layer thickness on shock structure in jets of 
arbitrary geometry - A finite difference technique has been developed to study 
the shock structure and hydrodynamic stability of jets with arbitrary exit geometry. 
A body-fitted coordinate system is used which is particularly suited for problems 
associated with ; bitrary exit geometry. This finite difference scheme includes the 
effect of a finite mixing layer thickness by using a realistic and continuous mean 
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Case 1: The image of the circle of R = n 


N 


6 

X 

y 

Error 

64 


0 

.749165 

.000000 

.002230 

128 


0 

.749791 

.000000 

.000590 

256 


0 

.749947 

.000000 

.00 139 

512 


0 

.749986 

.000000 

.000035 

64 


7r/4 

.528549 

.394927 

.004254 

128 


7r/4 

.529897 

.395195 

.001042 

256 


tt/4 

.530222 

.395262 

.000260 

512 


7t/4 

.5301 J2 

.395278 

.000007 

64 


7t/2 

.000000 

.562340 

.011924 

128 


tt/2 

.000000 

.559885 

.003110 

256 


7t/2 

.000000 

.559236 

.000787 

512 


7t/2 

.000000 

.559072 

.000197 


Case 2: 

The image of the circle R = 1 



N 


e 

X 

y 

Error 

64 


0 

.999086 

.000000 

.001826 

128 


0 

.999771 

.000000 

.000457 

256 


0 

.999943 

.000000 

.000114 

512 


0 

.999985 

.000000 

.000028 

64 


7t/4 

.706094 

.612594 

.001068 

128 


7t/4 

.706856 

.612424 

.000270 

256 


7t/4 

.707044 

.612385 

.000068 

512 


tt/4 

.707091 

.612375 

.000017 

64 


7t/2 

.000000 

.867770 

.004034 

128 


7t/2 

.000000 

.866468 

.001023 

256 


7t/2 

.000000 

.866136 

.000256 

512 


7t/2 

.000000 

.866053 

.000064 


Case 3: 

The value of 




N 






64 


.701323 




128 


.701456 




256 


.701489 




512 


.702497 






COMPUTATIONAL DOMAIN 



£>—« 


I 


PHYSICAL DOMAIN 


o 



>- 


profile. The method and the results for circular and elliptic jets are described 
below. 


Figure 4 shows the regions of a typical jet cross-section and the curvilinear 
coordinate system. The potential core region is bounded by the curve C i and an 
(s,n) coordinate system is used to represent it. The mixing layer is bounded by the 
curve Ci. The coordinates s and n are measured along and normal to the curve 
Ci respectively. The radius of curvature, R of curve Ci is a function of s only. 
Thus, the equation describing the curve C\ is all that is required to describe the 
coordinate system. 


The governing equations for a compressible, inviscid fluid including the curva- 
ture effects are given by, 


| + (1 + n/R))-' fj (pu) + + (1 


+ | -(pw) = 0 . 


du pu du du du puv 1 dp 

P dt "*"(! + n/R)) ds P dn PW dz (l + n/R)) R (1 + n/R)) ds 


dv pu dv dv dv pu 2 dp 

P dt (1 + n/R)) ds ^ pV dn + pW dz (1 + n/R)) R dn 

dw pu dw dw dw dp 

P dt + (1 + n/R)) ds + Pl dn + pW dz dz' 


pC v 


dT u dT dT dT_ 

dt (1 + n/R)) ds dn W dz 


dp u dp dp_ dp] 

dt (1 + n/R))ds V dn W dz j" 


Where (u,v,u>) are the velocity components in the (s,n,z) directions respectively, 


and p and T are the thermodynamic pressure and temperat r e respectively. 
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The general form of separable solution for the fluctuating pressure, when the 
mean velocity and density are taken to be independent, locally, of the axial distance 
is given by: 

p(s, n, z, t ) = p(s, n) exp[f [kz — w£)], 
where k is a wavenumber and oj is the radian frequency. 


The pressure fluctuation then satisfies the Rayleigh equation, 


dP 


+ 


dp 


dP 


dn 2 (1 + n/R))dn (l + n /R)) 2 ds 2 


2k 

n 


dU dp dU dp 


- (k 2 - pMfn 2 )p = 0, 


(1 + n/R)) 2 ds ds ' dn dn 
where, fi = (u>—kU) and the flow is assumed to be isothermal. This latter restriction 
may be readily removed. 


For an arbitrarily shaped jet U is a function of both s and n. Thus a separable 
solution in s and n is generally not available and one resorts to a numerical solution 
in the mixing layer. Along the interior edge of this region and in the potential core 
of the jet, the mean velocity is a constant. In the ambient fluid surrounding the 
jet the properties are also uniform. Thus on the inner and the outer boundaries 
of the mixing layer, a separable form of solution may be obtained. In general, 
a separable form of solution cannot be obtained in the ( s,n ) coordinate system. 
Thus, the general separable solution is obtained in (r, 0) polar coordinates and 
the relationship between the (r, 6) and («s,n) coordinates is used to transform the 
solution to the (s,n) coordinate system. This doesn’t pose any extra problems while 
dealing with arbitrary exit geometries. 

2.3.2 Numerical Calculations - The numerical method used here is identical to 
that used earlier, ref. 5, and hence, the derails of the method are not given here. The 
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only difference is that we are now integrating along lines normal to the potential 
core and not along radial lines. 

The mean velocity profile is taken to be, 

U(n ) = Uj exp[— In 2 (77 ) 2 ] , where 77 = n/b. 

b is the local half-width of the jet mixing layer. The relationship between h , the 
potential core radius, and b is obtained from an integral form of the axial momentum 
equation. 

Figures 5 , 6 and 7 show the variation of the axial wave number (lowest), which 
is inversely related to the shock spacing, as a function of mixing layer thickness and 
Mach number. Figure 5 shows the result for a circular jet at Mj = 1.4 and it can be 
seen that it compares well with the earlier work. Figure 6 compares the results for 
an elliptic jet of aspect ratio = 2 . For small values of 6 the spacing approaches that 
given by the vortex sheet model. Here, comparison is made with the results from 
the vortex sheet model and for b = 0.01 and 0.001. The agreement is very good 
for the Mach number range considered. Figure 7 shows the results for a circular 
jet and an elliptic jet of aspect ratio = 2 . In all these cases, only the lowest wave 
number associated with the axisymmetric mode has been sought. 

The results obtained here using the ( s,n ) coordinate system are in good agree- 
ment with the earlier work. The advantage of this method is that it is particularly 
suited for problems posed by arbitrary exit geometry. It should be noted that an 
inviscid model is used here and thus the calculations do not reflect the dissipative 
effect of the shear layer turbulence. The next stage of analysis is to include this 
viscous effect which would simply convert the equation satisfied by the pressure 
fluctuations from the Rayleigh equation to an Orr-Sommerfield equation. 
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Figure 5: Variation of Shock Spacing in Circular Jet with Mixing Layer Width 
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Figure 6: Variation of Shock Spacing in Elliptic Jet with Mach Number. 
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Figure 7: Variation of Shock Spacing in Elliptic and Circular Jets with Mixing 
Layer Thickness. 



The method described here can be used for other jet geometries, but at present, 
due to lack of experimental mean flow data at high speeds, it has not been imple- 
mented for, say rectangular and triangular jets. 

3. Turbulence Closure Scheme 

This part of the work is concerned with the prediction of the mean flow prop- 
erties of turbulent circular jets. A wave model is used to simulate the large-scale, 
coherent structures that dominate not only jets but also other free shear flows like 
mixing layers and wakes. Once the modeling procedure has been developed it will 
be extended to more complex geometries using the methods described in the pre- 
ceding sections. At the present stage we are focussing on the validation of the wave 
model and the development of an in-depth understanding of the physical character- 
istics of the wave structures. This is being assisted by studying a two-dimensional 
incompressible free mixing layer. This will provide guidelines for more complex ap- 
plications of the model. The equations used in the present wave model formulation 
are the continuity and the momentum equations of the long time-averaged proper- 
ties for the mean flow, Z7 t -, and of the wave-like large-scale structures, U{. These 
are 


dUi dui 

dx{ dxi 


Uj 


dUi 


+ 


dxj dxj 


■U{Uj + 


dxj 




where 


1 dP ' d 2 Uj 

g dxi dxjdxj 



Ti 


oo. 
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The nonlinear wave development may be simulated by introducing an ampli- 
tude function A(x) and describing the large-scale coherent structures in wave-like 
form: 

{u,u} = A(x) 9?{[<^(j/),^(j/)] exp[i{ax- u><)]} . 

The amplitude of the wave-like disturbances is determined by the kinetic energy 
equation for the large-scale motions, 


U< 


dk 

dxj 


= —UiUj 


_dUj 

dxj 


- (- < >) 


du{ 

dx. 


— — — (u t - < u'-u'- >) + viscous terms 
dxj 1 > 


where 

1 f T » 

<q>=7fr q dt, T 2 «Ti 
■t 2 Jo 

Our previous calculations, ref. 6 have shown that the influence of the small-scale 
structures cannot be neglected. This is especially important on the low-velocity 
side of the mixing layer. It has also been shown experimentally ref. 7 that the 
contribution to the long time-averged Reynolds stresses from the large and small 
scale components may be equally important in fully-developed mixing layers. A 
split-spectrum approach was thus adopted ref. 8, which assumes that some of the 
turbulent kinetic energy is contained in the small-scale motions. A computer pro- 
gram implementing an interactive approach that links the long time-averaged flow, 
the wave-like, large-scale motions, and the small-scale fluctuations has been de- 
veloped. The small-scale, non-coherent stresses may be modeled initially using an 
eddy-viscosity form, lation, that is, 
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The starting mean velocity profile was taken from experimental tesults for a 
free mixing layer. It was found that in the long tai region of the velocity profile 
the negative production of the large-scale structures prevented the flow from ac- 
celerating and stopped the mixing region from growing. Figure 8 gives the mean 
velocity profiles in the streamwise direction. The deceleration effect can be clearly 
seen. Similar results were found for a two-stream mixing layer. For a free mixing 
layer, the acccuracy of the measured mean flow data on the low-speed are subject 
to error due to the rapid variations in the instantaneous flow direction. Thus it is 
possible that the waves driven by the mean flow see a “mean” different from that 
obtained from a long time-average, especially at places where the fluctuations may 
actually modify the “mean” flow as they pass by. 

We also found that the amplitude of the waves grew exponentially. It is known 
from experiments, however, that the amplitude grows exponentially for only a short 
distance downstream before saturating. In our models it is the residual stress tensor 
7t/ which appears in the wave kinetic energy equation in the form 


d u, 



that is responsible for the draining of the energy from the waves. The residual 
stresses were simulated by the eddy-viscosity model. The constant in the model 
was taken to be the same as that in the eddy-vicosity model for the small-scale 
Reynolds stresses. It was found that the term thus calculated was far too small to 
represent the amount of energy that would be transfered from the large-scale to 




Figure 8: Variation of Mean Velocity with Axial Distance. 



that is responsible for the draining of the energy from the waves. The residual stresses 
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where C is a constant and l is a local length scale. This is based on dimensional reasoning 
and the assumption that the 7 t y scale with the kinetic energy, k, of the wave motions. Part 
of the turbulent kinetic energy is then taken away from the motion of the wavy compo- 
nents and limits their amplitude. Figure 9 shows the development of the mean stream-wise 
velocity component calculated using this model. We also compared the variations of the 
wave amplitude in these two cases. The results are shown in figure 10. The amplitude 
function initially grows exponentially and then flattened. This agrees with experimental 
findings, ref. 8. Thus, a good model for these residual stresses is crucial since the coher- 
ent structures dominate the flow development and their amplitude must be determined. 
Unfortunately, little experimental data is available regarding these stresses. 


Several approaches are now being developed to model this energy transfer mechanism. 
On-going efforts also include using the kinetic energy equation for the small-scale motion 
to estimate the velocity scale in the eddy viscosity model. The equation for the small-scale 
turbulent kinetic energy is 
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/c 3 / 2 



where C is a constant and l is a local length scale. This is based on dimensional 
reasoning and the assumption that the scale with the kinetic energy, k, of the 
wave motions. Part of the turbulent kinetic energy is then taken away from the 
motion of the wavy components and limits their amplitude. Figure 9 shows the de- 
velopment of the mean stream-wise velocity component calculated using this model. 
We also compared the variations of the wave amplitude in these two cases. The re- 
sults are shown in figure 10. The amplitude function initially grows exponentially 
and then flattened. This agrees with experimental findings, ref. 8. Thus, a good 
model for these residual stresses is crucial since the coherent structures dominate 
the flow development and their amplitude must be determined. Unfortunately, little 
experimental data is available regarding these stresses. 

Several approaches are now being developed to model this energy transfer 
mechanism. On-going efforts also include using the kinetic energy equation for 
the small-scale motion to estimate the velocity scale in the eddy viscosity model. 
The equation for the small-scale turbulent kinetic energy is 
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Note that the equation must be solved simultaneously with the mean flow 
equations. Initially the modeling of each term in the equation is performed in a 
similar way to that of traditional closure schemes. 
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